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Quantum walks and quantum search on graphene lattices 
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Quantum walks have been very useful in developing search algorithms in quantum information, 
in particular for devising of spatial search algorithms. However, the construction of continuous-time 
quantum search algorithms in two-dimensional lattices has proved difficult, requiring additional 
degrees of freedom. Here, we demonstrate that continuous-time quantum walk search is possible in 
two-dimensions by changing the search topology to a graphene lattice, utilising the Dirac point in 
the energy spectrum. This is made possible by making a change to standard methods of marking a 
particular site in the lattice. Various ways of marking a site are shown to result in successful search 
protocols. We further establish that the search can be adapted to transfer probability amplitude 
across the lattice between specific lattice sites thus establishing a line of communication between 
these sites. 

PACS numbers: 03.67.Hk, 72.80.Vp, 03.65.Sq, 03.67.Lx 


I. INTRODUCTION 

Quantum walks have been a particularly fruitful field 
of research in quantum information going back to ideas 
from Feynman [I], Meyer [2], and Aharonov fi . One 
focus of recent work on quantum walks is their applica¬ 
tion to spatial search algorithms. It is well-known that 
Grover’s algorithm H El for searching on unstructured 
databases offers a quadratic improvement in search time 
over classical models. Grover’s algorithm is not designed 
to search physical systems where only local operations 
are possible and quantum walk algorithms have been em¬ 
ployed for this purpose. In developing such algorithms, 
the major consideration is the search time attempting to 
match the quadratic improvement over the classical case 
offered by Grover’s algorithm. The first spatial search al¬ 
gorithm to do this, using the standard model of quantum 
walks, was developed by Shenvi, Kernpe and Whaley [5] 
for searches on a hypercube. 

While there exist discrete-time searches on d- 
dimensional cubic lattices which are faster than classical 
searches for d > 2 (7] , effective continuous-time quantum 
searches only exist for d > 4 jBj or else they require ad¬ 
ditional memory (in the form of spin degrees of freedom) 
in order to improve their search time in lower dimen¬ 
sions m- In [TO] , we demonstrated that effective searches 
over two-dimensional lattices may be achieved in an ar¬ 
guably simpler way which does not require extra degrees 
of freedom, and could, therefore, be viewed as more effi¬ 
cient. This is achieved through the choice of a different 
lattice, specifically, a honeycomb lattice which is the un¬ 
derlying lattice structure of carbon atoms in the material 
graphene. 

The association with graphene is important as, al¬ 
though we first study a purely theoretical problem in 
quantum information, the use of a graphene lattice also 
offers a potential physical realisation. We thus envis¬ 
age using the quantum walk and quantum search algo¬ 
rithm framework to investigate the effect of perturbations 
on the dynamics on graphene and other carbon struc¬ 
tures. This offers not only the possibility for demonstrat¬ 


ing two-dimensional continuous-time quantum searching, 
but also paves the way for looking for novel effects in the 
material graphene. 

In this paper, we offer a detailed account of the theory 
and numerical results thus expanding on the findings as 
presented in |lf)| . Generalisations of the results in [10] 
have been given in m, where a approach to solving the 
problem has been taken through encoding extra degrees 
of freedom into crystal lattices. Recent experiments on 
artificial microwave graphene [15] have shown that ad¬ 
ditional site perturbations, as discussed in Sec. m can 
be used to create a search protocol. Searching on hon¬ 
eycomb lattices in discrete-time setting has been consid¬ 
ered in [13] ; however, no improvement over discrete time 
searches on cubic lattices was found. 

The paper is structured as follows: In Section [TT] we 
give an introduction to the formalism of continuous-time 
quantum walks and the construction of quantum walk 
search algorithms. We also explain why previous search 
algorithms struggled on lower-dimensional lattices and 
why graphene offers a solution. Section m will detail 
the relevant properties of graphene and set-up the no¬ 
tation we shall use throughout. Sections ||V] & [V] con¬ 
tain our main analytical results, where we detail the 
specifics of our search algorithm and offer an analysis 
of the search running time and success probability. We 
shall then show in Section [VT| how this search protocol can 
be adapted to demonstrate novel communication setups. 
Sections IVIII & IVIIII contain numerical work demonstrat¬ 
ing the possibility of using alternative methods of mark¬ 
ing to create search behaviour and other carbon nanos¬ 
tructures. We conclude with a review and discussion of 
our results in Section m 


II. QUANTUM WALKS AND SEARCHING 

Continuous-time quantum walks (CTQW), first de¬ 
fined in [14], are the quantum analogue of continuous¬ 
time Markov chains. They are defined purely on the 
state space, that is, the Hilbert space spanned by 
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the states | j) which represent the j th site of the lattice. 
Thus, the time-evolution of such systems is defined by 
the Schodinger equation 


d 

dt 


CXj (t) 


N 

(*) > 

1=1 


(i) 


where aj = (j \ ip (f)) is the probability amplitude at 
the j th vertex of a system described by the state vector 
| if) (i)), and H is the Hamiltonian describing the connec¬ 
tivity of the lattice. Note that we are using a dimension¬ 
less description setting e.g. h = 1. 

The dynamics of a quantum walk over a network is 
defined by the nature of the interaction between con¬ 
nected sites. Therefore, the Hamiltonian is generally 
constructed from the adjacency matrix of the underly¬ 
ing lattice. The adjacency matrix A is defined as 

_ f 1 if j and l are connected 

_ \ 0 if j and l are not connected. ' ' 

Typically, the Hamiltonian is chosen as H = e^I + vA. 
The parameter v determines the coupling strength be¬ 
tween connected sites and the parameter cd is an on-site 
energy that only enters the dynamics in a trivial way and 
thus can be set to a desired value. If v = — 1 and en is 
equal to the valency of the lattice H is the discrete Lapla- 
cian. This form of Hamiltonian is closely related to the 
tight-binding model for condensed matter systems US¬ 
As first explained in puna, a quantum walk is trans¬ 
formed into a search protocol by introducing a localised 
perturber state, forming an avoided crossing in the spec¬ 
trum of the search Hamiltonian between an unperturbed 
eigenstate and the localised perturber. Thus, initialising 
the system in the unperturbed eigenstate involved in the 
crossing and allowing the system to evolve in time, one 
finds that the system rotates into the localised perturber 
state. 

The first CTQW search over d-dimensional cubic lat¬ 
tices jH] introduced the localised perturber state as a pro¬ 
jector onto a single site w, resulting in the search Hamil¬ 
tonian 


H t = —yA + |u>) (u>| . (3) 

Here, 7 is a parameter governing the strength of inter¬ 
actions between sites in the lattice. This parameter is 
chosen carefully 0 such that the perturber state \w) is 
brought into resonance with the ground state of the un¬ 
perturbed Hamiltonian, the uniform superposition |s). 

Let us assume for the moment that all other unper¬ 
turbed eigenstates are energetically sufficiently separated 
from the ground state - we will come back to this assump¬ 
tion later. Then there will be an avoided crossing of two 
eigenvalues corresponding to the perturber state and the 
ground state. Perturbation theory estimates that the en¬ 
ergy splitting of these resonant states is proportional to 
the overlap of the localised perturber state |w) and the 
uniform ground state |s), which scales as 

AE ~ |(ic|s)| ~ A -1 / 2 , (4) 


and that the corresponding eigenstates are of the form 
(|s) ± \w))/\/2. By preparing the system initially in |s) 
and allowing it to evolve for a set period of time T = 
tt/AE oc y/N , a measurement of the system will result in 
the state |w) being measured with high probability. This 
explains the speed-up of the quantum walk search. The 
Grover algorithm and its speed-up can be understood in 
analogous terms and one can show that all other states 
are indeed energetically very well separated. 

In the present case this separation only holds above 
a critical dimension d c = 4. A simple argument for the 
critical dimension can be obtained by comparing the scal¬ 
ing behaviour of the energy splitting AE ~ A -1 / 2 with 
the energy separation of the first excited state from the 
ground state in the unperturbed lattice. For a cubic lat¬ 
tice we have a quadratic dispersion relation which allows 
us to estimate the energy separation 

E{k) - E 0 ~ |fc | 2 ~ N~ 2/d 

(as |fc| ~ A _1 / d for the first excited state in a d- 
dimensional lattice). For d > 4 one then has A E/(E(k) — 
Eq) —> 0 as N —> 00 and a detailed analysis indeed proves 
that the quantum walk search works with optimal speed¬ 
up 0 - At the critical dimension d = 4 the two energy 
scales scale in the same way - the detailed analysis shows 
that the search still works but the dynamics is more com¬ 
plicated due to the interference of excited states. While 
there is still a speed-up it is only almost optimal; the 
optimal search time T oc A 1 / 2 gets multiplied with a 
logarithmically increasing factor. For the experimentally 
relevant regimes of 2- and 3-dimensional cubic lattices, all 
states participate in the dynamics as any avoided crossing 
gets dissolved completely as the number of sites grows. 
As a result no speed-up over classical searches is found. 

The above estimate offers a simple way how to reduce 
the critical dimension. If one can construct a search 
around a uniformly distributed state at an energy Eq 
where the dispersion relation is conic (linear in |fc|), then 
in d dimensions 


E — Eq ~ \k\ ~ N 1/d 

which results in a critical dimension d c = 2. In (93 
this was implemented using a modified Dirac Hamilto¬ 
nian. However, this requires the addition of a spin de¬ 
gree of freedom, essentially a doubling of memory, with 
the added complication that it is not immediately clear 
how such a system would be physically realised. 

Instead, our solution here is to change the lattice from 
a square (cubic) to a graphene lattice. This change of 
topology automatically implies the first step in our so¬ 
lution, as one of the important electronic features of 
graphene is the conic dispersion relation around the Dirac 
energy. This arises naturally from the tight-binding de¬ 
scriptions of graphene. Note that graphene has the crit¬ 
ical dimension d = 2. Indeed our construction as pre¬ 
sented in m has an almost optimal speed-up with log¬ 
arithmic corrections that need to be evaluated in a de- 
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FIG. 1. Left: Graphene with lattice vectors Oj , 2 , translation 
vectors and unit cell (dashed lines). Right: Reciprocal 
lattice with basis vectors , 2l symmetry points T, K _, K ', M 
and first Brillouin zone (hexagon). 


tailed analysis that goes beyond the simple perturbative 
description given above. FIG. 2. (Color online) Dispersion relation for infinite 

graphene sheet (en = 0 and v = —1). 



III. RELEVANT PROPERTIES OF GRAPHENE 

Graphene is a single layer of carbon atoms arranged 
in a honeycomb lattice. The lattice is bipartite with two 
sublattices, labelled A and B, and a unit cell containing 
two carbon atoms. The spatial and reciprocal lattices 
are shown in Figure [l] The primitive vectors describing 
the lattice are a-^ 2 ), such that the position of a unit cell 
in the lattice is given by R(a,(3) = aa\ + (3a 2 - We use 
dimensionless units in space where the distance between 
nearest neighbor sites is a = 1. The reciprocal lattice 
shows two important points, the Dirac points K_ and K* 
at the two inequivalent corners of the Brillouin zone. 

The energy spectrum of electrons in graphene was first 
derived by Wallace [T8] when considering the band struc¬ 
ture of graphite using a tight-binding Hamiltonian 

H = e£)l + vA (5) 


with a reduced density of states, a necessary feature for 
the creation of the search dynamics. 

As the lattice possesses a translational symmetry the 
Hamiltonian can be solved using linear superpositions 
of Bloch functions over both sublattices. As a basis we 
use the orthonormal states {|a,/3) , |a,/3) } to denote 
states on either the A or B sublattice in the cell at posi¬ 
tion R ( a , (3). For the majority of what follows (except in 
Section |VIII ), we will focus on finite-sized lattices with 
assumed periodic boundary conditions along the axes of 
both primitive vectors so that the topology of our lattice 
is a torus; that is, our wavefunction is of the general form 


E ra 

a=l 2.^8=1 


•ip\ 


Our 


AP) A +il>^p\a,/3) B y 

boundary conditions imply that the state vector must 

satisfy = V’a+m./S = ^/ 3 +n where m , n denote the 

period of the lattice. Thus, the wavefunctions on the 
torus take the Bloch function form 


where tp is the on-site energy (which we will identify 
as the energy of the Dirac points) and v the hopping 
strength (both dimensionless in our setting). The tight- 
binding model for graphene and the derivation of the 
solution are well-known PM and give rise to the dis¬ 
persion relation 


e (k) = e_o± 


( 6 ) 


V 


N 


1 + 4 cos 2 




shown in Figure [2] for an infinite graphene lattice. 

As there are two atoms per unit cell the spectrum has 
two branches, the upper branch being the conduction 
band and the lower the valence band, which meet at the 
corners of the Brillouin zone, the A-points. The energy 
at the AT-points is £d which we name the Dirac energy. 
It is around these points, that the behaviour of the spec¬ 
trum is conical, that is, 


e{k) 


. V3 
£D±V T 



, V5 
eD±v T 


\Sk\ , 


(7) 


= £ 

(“. 0 ) 


A_e ik -^\a,P) A 


c ik-(R( a ,P)+6 ,) | a;/3 jB 


( 8 ) 


where N is the number of sites in the lattice, k is the 
momentum, and C ( k ) is a relative phase contribution 
dependent upon whether the state belongs to the conduc¬ 
tion or valence bands; it can be calculated by explicitly 
working through the tight-binding model. 

Application of the periodic boundary conditions results 
in the following quantised momenta 


7,: — 


2np 1 f Airq 

m V3 V n 


(9) 


where p £ {0,1,... m — 1} and q £ {0,1,... n — 1}. In 
the following and whenever we consider quantum walk 
dynamics on a torus, the number of cells in each direction 


is generally chosen to be the same, that is, m = n = y T- 
This choice is purely for simplifying the notation; alter¬ 
native torus dimensions are possible as are other choices 
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of boundary conditions corresponding to alternative car¬ 
bon structures (e.g. nanotubes or a graphene sheet) as 
will be demonstrated in Section IVIIII For our choice of 
torus dimensions, we find there are momenta equal to 
the IY-points and, consequently, eigenstates with ener¬ 
gies equal to the Dirac energy when both m and n are 
some multiple of 3. 

In fact, using the quantised momenta in Eq. ([9]) ob¬ 
tained for periodic boundary conditions, we find that 
there are four degenerate eigenstates with an energy that 
coincides exactly with the Dirac energy when m and n 
are both multiples of 3. These four states, known as the 
Dirac states, can be constructed to live only on one of 
the sublattices, and are given by 



V) A(B) =\fl £ e' w ‘ +f) K/A (B> ■ (io) 

V («,/3) 

where a = 0 for states on the A-sublattice or a = 1 for 
states on the H-sublattice. 


IV. QUANTUM SEARCH ON GRAPHENE - 
REDUCED MODEL 

In this section we will describe how a site is marked 
and will derive our optimal search starting state. Our 
approach will be to analyse the system’s spectrum and 
its dynamics in a reduced Hamiltonian model involving 
only the relevant states from the avoided crossing. In the 
next section we will then validate this approach with a 
more detailed analysis, using the results obtained here as 
an initial guide. 

As already established, we introduce the localised per- 
turber state in a region of the spectrum with a conic 
dispersion relation and thus a low density of states. Sim¬ 
ply altering the on-site energy in Eq. ^ as done in [ 8 ] 
does, however, not work here. As the on-site energy e£> 
and the Dirac energy are equal, one finds that the per¬ 
turbation only interacts with the Dirac states in the limit 
of zero perturbation strength, returning the unperturbed 
lattice. Therefore, we choose an alternative perturba¬ 
tion method: namely, we modify the coupling strength 
between the site we wish to mark and its neighboring 
vertices. We focus in this section on changing the cou¬ 
pling to all three neighboring vertices of a particular site 
equally. Our choice of perturbation matrix W to mark 
the A-type vertex ( a 0 ,P 0 ) A is then 

W = \/3 \t) (a 0 , p 0 \ A + Vs \a 0 , p a ) A (£\ , (11) 

where the state \£) is the symmetric superposition over 
the three neighbors of the perturbed site 

\£) = (jcG ~ 1, Po) + Vo — 1) Po + 1) + Vo, Po) ^ • 

(12) 
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FIG. 3. (Color online) Spectrum of H 7 in Equation ( |13| ) as 
a function of 7 for a 12 x 12 cell torus (N = 288). The 
spectrum is symmetric around to = 0. Inset: Scaling of the 
gap A = E+ — E _ (dots) and curves Ci/V~N (solid blue), 
C 2 /VN log N (dashed red) for comparison. 


This leaves us with the search Hamiltonian 

H 7 = - 7 A + W, (13) 

where 7 is a free parameter. In what follows, we always 
set the on-site energy en =0. Considering our search 
Hamiltonian, we can see that setting 7 = 1 corresponds 
to a coupling strength of v = 0 from the perturbed vertex 
and its nearest-neighbors; our perturbation essentially 
removes the marked vertex ( a 0 ,P 0 ) A from the lattice. 
Note that vacancies are a common, naturally occurring, 
defects in graphene lattices [ 22 ] , 

In order to establish the critical value of 7 , we numer¬ 
ically calculate the spectrum of H 7 as a function of 7 , 
plotted in Figure[3]for a torus of dimensions m = n = 12. 
As W is a rank-2 perturbation, we see in Figure [3] two 
perturber states working their way through the spectrum 
to an avoided crossing around eu = 0 when 7 = 1 , that 
is, when the perturbed vertex is removed from the un¬ 
derlying lattice. 

As well as establishing a critical value for 7 , we also 
find from this figure the states involved in the avoided 
crossing: the four Dirac states at the Dirac energy and 
the two perturber states which form our perturbation 
matrix W. However, with further consideration, we may 
reduce the number of states involved further. As we 
have removed the site (ao, Po) from the lattice it can no 
longer interact with the rest of the lattice and the cor¬ 
responding state drops out. Also, by direct calculation 
one finds that the U-type Dirac states do not interact 
with the perturbation, that is, W \K) B = W \pV) B = 0. 
Thus, they remain an eigenstate of the search Hamilto¬ 
nian and do not interact with the perturbation. We are 
left with three states taking part in the avoided crossing: 

{\k) a ,\k') a ,\m- 

We reduce our search Hamiltonian in this three state 
basis at the critical point 7 = 1 to obtain the following 
reduced Hamiltonian describing the local dynamics at the 
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FIG. 4. (Color online) Search on 12 x 12 cell graphene lattice 
with a triple-bond perturbation, using starting state |s) from 
Equation i m For tori with m = n the dynamics at each 
neighboring site is the same so only one is shown. 


avoided crossing 


H = 



0 

0 



(14) 


with (a 0 + 2/3 0 ) and = vf (2 a 0 + /3 0 ). 

This reduced Hamiltonian has eigenvalues E± = 


±2 



, Eq = 0 , and eigenvectors 


7±) = \ | K) a + e~ iv ° | K') A ± 72 |£)) (15) 

7o) = ^ (e-^° | K) a - eT iv ° \K') A ) . (16) 

Using the eigenvectors of the reduced Hamiltonian, we 
can construct a search starting state which is a superpo¬ 
sition of Dirac states 


17 


1 

72 


( 7 +) + 7 -)) 




(17) 


Allowing our search starting state |s) to evolve under 
the reduced Hamiltonian we find 


\4,(t))=e- Mt \a) 

1 / 


72 




— iE+t 


7 + 


—iE—t 


7 - 


(18) 


cos I s ) — *sin (^+i) K) , 


so that our system rotates from our Dirac superposition 
|s) to a state which is localised on the neighbors of the 


perturbed vertex |£) in a time t = f y y - Thus, we find 
a O I VN) search time, a polynomial improvement over 


the classical search time. 

We plot in Figure [4] a numerically calculated search 
for a lattice with N = 288 sites. The system has been 


prepared in |s) and allowed to evolve under the full search 
Hamiltonian from Eq. (13) at the critical value 7 = 1. 


One finds that the system localises on the neighbors of 
the marked vertex and, for this particular system size, 
the probability for a measurement of the system to return 
one of the neighbors is around 45% at its peak. This is 
orders of magnitude higher than the average probability 
to measure other vertices in the lattice, lOO/IV, which 
for this system size is around 0.35%. The localised state 
interacting with the spectrum under the full Hamiltonian 
has a tail into the rest of the lattice; this leads to a loss of 
probability to be found at the neighboring vertices below 
100 %. 

It is worth emphasising that the marked vertex plays 
no role in the search - it is removed from the lattice. 
Rather, we force the system to localise on the nearest- 
neighbors, essentially making the marked vertex conspic¬ 
uous by its absence. This feature differentiates our algo¬ 
rithm from previous searches, which have generally con¬ 
centrated on localising directly on the marked vertex. 
An approach considering the neighboring vertices was 
detailed in |27| , however, the focus was here still on lo¬ 
calising on a single site but with classical post-processing 
steps considering the tail of the localised state extending 
into the neighborhood around the marked vertex. 

At this point, there are a number of issues which need 
to be mentioned. The first is that our starting state in 
Eq.@ contains information about the marked vertex in 
the form of the relative phase between the Dirac states. 
However, this phase can only take three different values. 
Therefore, we have three possible optimal starting states 
for an A-type perturbation and the same number for B- 
type perturbations; there are thus in total six possible 
optimal states. As not all of these states are orthogonal, 
we only find an increase of necessary runs for a successful 
search by a factor of 4; this additional overhead is inde¬ 
pendent of N and, therefore, does not affect the overall 
time complexity of the search. The particular represen¬ 
tation of the Dirac states in Eq. (10) was chosen in such 


a way as to make the calculations and conceptualisation 
of the system dynamics easier. In reality, constructing 
starting states which exist on one sublattice only, experi¬ 
mentally is extremely difficulty; rather, in an experiment, 
one is likely to excite a superposition of Dirac states, re¬ 
ducing the success probability, on average, by a factor of 
1/4. 

Our search is based on the conic dispersion relation in 
the spectrum and the O (l/VN^ scaling between succes¬ 
sive energy levels in the linear regime, giving rise to the 
O ( VN ) search time found in our reduced model. How¬ 


ever, previous searches in continuous-time using a mod¬ 
ified Dirac Hamiltonian [9] or in discrete-time 7] have 
found logarithmic corrections to the search time. As we 
have seen in Eq. (18) the search time is related to the 


energy gap at the avoided crossing. In the inset of Fig¬ 
ure [3] we have numerically calculated the scaling of the 
gap in the spectrum of the full search Hamiltonian and 
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we indeed find evidence for a logarithmic correction. The 
lack of a logarithmic term in the search time for our re¬ 
duced model is due to the neglecting contributions from 
the rest of the spectrum. 

In what follows we establish a more accurate estimate 
of the running time and success probability by working 
with the full Hamiltonian. We also justify the findings 
of our reduced model: while this model is insufficient to 
estimate the finer details, such as logarithmic correction 
terms, the accurate calculations show that our search al¬ 
gorithm indeed takes place mainly in a two-dimensional 
subspace spanned by the Dirac states and a state lo¬ 
calised on the neighbors of the marked vertex. 


V. QUANTUM SEARCH ON GRAPHENE - 
DETAILED ANALYSIS 

We follow here closely our derivation as presented in 
m which builds on ideas given in |5] for a similar a 
similar calculation for regular rectangular lattices. We 
focus first on the search success amplitude, that is, 


(t\e 


-iUT 


| start) = ^2 (^1 V’a) (VO start) e 
I W> 


—iE„T 


(19) 


where T is the search time, that is, the time our search 
probability reaches a maximum, and our search Hamil¬ 
tonian H is given by Eq. (131 with 7 = 1, that is, 


H = —A + V3 \£) (a 0 , Po\ + VH |a„, /3 0 ) {£\ ; (20) 

here, | ip a ), E a are the eigenstates and eigenenergies of H. 
We assume, without loss of generality, a specific starting 
state | start) where the marked vertex is chosen such that 

e >f (a 0 +2^ 0 ) _ i i eac ii n g to 


I start) = -j= (| K) + | K')) . 


( 21 ) 


In the following, we suppress the sublattice superscript 
as the analysis is the same regardless on which sublattice 
the perturbation lives. We also denote e (k) the positive 
eigenenergies of —A from Eq. § and set cd = 0. 

For eigenstates \ip a ) with eigenenergies E a that are not 
in the unperturbed graphene spectrum (E a ^ e(k) for all 
points k in the dual lattice), we may rewrite Eq. (201 in 
the form 


m = ^Ra(E a + A) 1 \a 0 ,p 0 ) 


( 22 ) 


of H which are at the same time eigenstates of — A with 
the same energy. 

This can be seen in the following way. We first consider 
an unperturbed eigenstate | ip°) such that —A \ip°) = 
E a \ip°). Let us assume that there is an eigenvector 
\ip a ) of the search Hamiltonian with the same eigenen- 
ergy E a , that is, H| i(> a ) = E a \ip a ). Considering the 
matrix element (ip°\ H | ip a ), we find (ip°\ £) (a Q , /3 a \ ip a ) + 
(ip° I ct 0 ,/3 0 ) (£\ ip a ) = 0. As |a 0 ,/3 0 ) is an eigenvector of 
the search Hamiltonian we know that (a c ,/ 3 0 \ ip a) = 0. 
This leaves us with {ip°\ a 0 ,/3 0 ) {£\ ip a ) = 0. As the un¬ 
perturbed eigenstate | ip°) is simply a Bloch state we know 
{ip° | a 0 ,/3 0 ) 7^= 0. Thus, we obtain {£\ ip a ) = \£Ra = 0. It 
is then clear that eigenstates of the search Hamiltonian 
whose eigenenergies remain in the spectrum of —A do 
not play a role in the time-evolution of the search. 

We now derive an eigenvalue condition for those per¬ 
turbed eigenvalues which are in the spectrum of H. Using 
the orthogonality of eigenstates (a 0 , j3 a \ ip a ) = 0, Eq. (22) 
leads to 

y/SRa ( a o, Po\ (E a + A) -1 \a 0 , f3 0 ) = 0. (23) 

By expressing (a 0 ,/3 0 \ in terms of the eigenstates of the 
unperturbed walk Hamiltonian —A, we may write this 
as a quantization condition 


with 




F (E a ) = 0 


+ 


E — e(k) E + e(k) 


(24) 


(25) 


were N is the total number of sites in the lattice. Eq. (25) 


is written as two summations to incorporate the fact that 
the spectrum of —A as well as H is symmetric around 
E = 0. 


Choosing \ip a ) to be normalised (ip a \ ipa) = 1, Eq. (22) 
also implies 

3 Ra (c*o, Po \ (E a + A)- 2 \a„, (3 0 ) = 1, (26) 

which allows R a to be rewritten as 

1 


Ra = 


V3| F’{E a )\ 


(27) 


where = {£\ ip a ) and the phase of \ip a ) is chosen such 
that (£| ipa) > 0. 

At this point, we can remove several states from the 
summation in Eq. (19). Note that the basis state associ¬ 
ated with the marked vertex, \a 0 ,/3 0 ), is itself an eigen¬ 
state of the search Hamiltonian with H \a Q , /3 a ) = 0, and 
also, {£\ a 0 ,/3 0 ) = 0 so that \a 0 ,/3 0 ) does not contribute 
to the time-evolution. We may also remove eigenstates of —A. 


We may now rewrite the amplitude in Eq. (19) in the 
form 

e -iE a T 

{£\ e~ lH1 \start) = (a 0 ,/3 0 | start) ^ , , , , 

_ ( 28 ) 

where we have used the adjoint of Eq. (22). Note again 


that eigenstates of H, which are also in the spectrum 
of —A, do not contribute to the time-evolution of the 
search; it may be seen from the definition of F ( E ) in 
Eq. (25) that |F' (E a ) \ -A oo, where E a is in the spectrum 
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As we saw in the previous section, the avoided cross¬ 
ing is formed by two perturbers approaching symmetri¬ 
cally either side of the Dirac states with energy e D = 0. 
Thus, we concentrate on evaluating the contribution to 
the time-evolution from these perturbed eigenstates of H 
either side of the Dirac point and we label these states 
\ip±). In order to evaluate these contributions we cal¬ 
culate the perturbed eigenenergy E + and also derive a 
leading-order expression for F' ( E + ) (since the spectrum 
is symmetric, we have E + = —E- > 0). 

Using the definition of F ( E ) in Eq. (|25j), we esti¬ 


mate F (E + ) by separating out the Dirac points, where 
e ( K_) = e (lU) = 0, from the summations and then Tay¬ 
lor expanding the remaining terms at E = 0, to find 




(29) 


where the sums I n are given by 

1 1 


j VI Y 
ln ~ N 


hc*)] b J 


(30) 


As the unperturbed spectrum is symmetric only those I n 
with even n are non-zero; thus from now on we focus only 
on 1 2 k where k > 1. 

We stated earlier that the spectrum of graphene is well- 
approximated by a conic dispersion relation around the 
Dirac points. Thus, it is from around these points that 
the major contributions to the l 2 k sums arise. By ap¬ 
plying the linear approximation from Eq. (J7| and the 
momenta quantum numbers from Eq. ([9]) we may ap¬ 
proximate the 12 k summations as 


hk = 4v / 3Af fc_1 


+ 0 ( 1 ). 


E Z 2 {Sk, 2)+ ^ Z 2 {Sfc. 2) 

(p,q)eL (p,q)£L' 

(31) 


Convergence of the Epstein zeta function is well-known 
for k > 2 j2TJ , and leads to the bounds 

= 4v/3 (Z 2 (Sk, k) + Z 2 (S K ',k )) (34) 

for k > 2. A sharp estimate for 


I 2 = O (ln AT) . 


(35) 


is given in Appendix [A] In order to calculate an estimate 
for E_ |_, we truncate the Taylor expansion of F (E) in 
Eq. ( f29j ) at the I 2 term (the first term in the summa¬ 
tion), and apply the eigenvalue condition F (E + ) = 0. 
That is, we solve — I 2 E + = 0 and obtain the ap¬ 
proximation E+ ~ This solution also leads to the 

estimate F' ( E + ) ss —2 1 2 . 

We consider next whether our solution for E + lies in¬ 
side the radius of convergence of the Taylor expansion. 
We note that each term in Z 2 (Sx,k ) is smaller than 
the corresponding term in Z 2 ( Sk _, 2) for k > 2, and so 
it follows that ^(.SV, k) < Z 2 (Sk^, 2) for k > 2. This 
property of the Epstein zeta function and the estimate 
from Eq. (34) imply 


OO 

J2 ^ nE +”” 1 < 

n =2 


c 


NE 


Y2( NE ir 


n=2 


i.e. the infinite sum in Eq. (29) converges for E + < 


1 /\fN. Thus, for large N, our solution lies within the 
radius of convergence of our Taylor expansion in Eq. (291. 


In order to show that the dominant contributions to 
the search come from \i/j±), we need to establish the lead¬ 
ing order error term. All the 12k sums are positive so 
that, given the sign of all the terms in the Taylor ex¬ 
pansion in Eq. (291, the true value of E + > 0 has to be 
smaller than the estimate we have obtained. Thus, we 
write the true value of E + as 


Here Z 2 {S,x) is the Epstein zeta-function |2TJ 

Z 2 (S,x) = Y {SnP 2 + 2S V2 pq + S 22 q 2 )~ X , 

(p,q)e z 2 \(o,o) 

,( 32 ) 

for a real positive definite real symmetric 2x2 matrix S. 
For our purposes we use 

Sk = S K > = 4tt 2 , (33) 


which describes the spectrum close to the Dirac points. 
The linear approximation around both Dirac points K 
and K' is the same and, thus, the matrices Sk and Sk 1 
are equal. 

Our summation in Eq. (311 is over the rectangular re¬ 
gions L and L' of the lattice Z 2 . Both are centered on 
(0, 0) and have side lengths proportional to VN, however 
the center, (0,0), corresponding to the relevant Dirac 
point, is omitted from the summation. 


E l = mr A>Q ’ <36) 

with A > 0. One may rewrite F(E + ) = 0, using the true 
value of E + , to give 


I 2 A = 51 hnE 2 + n . (37) 

n —2 


We follow the same arguments as used for the calculation 
of the radius of convergence to obtain an upper bound 
for the summation. This together with the already es¬ 
tablished fact that E + is inside this radius for sufficiently 
large N, we get the following inequality: 


0 < NI 2 A < Cj2(NEl) n 

n =2 


CN 2 E\ 

l-NE'l 


o(/ 2 - 2 ). 




















So A = O^A" 1 ) = O ((In A)" 3 A" 1 ). Thus, we ob¬ 
tain 


^ = ^(l + 0((lnA)- 2 )) . 


(38) 


We use the same unperturbed walk Hamiltonian as be¬ 
fore, that is, Ho = — A and the same type of perturbation 
matrix as in Eq. (11), but now at two different sites in 
the lattice. Explicitly, our communication Hamiltonian 
is 


It also follows that 


H = -A + W s + W t , (45) 


F' (E±) =-2I 2 + 0 


In A 


(39) 


This shows that our perturbed eigenstates | if>±) have 
an O (l)-overlap with the starting state and are, there¬ 
fore, the relevant states to be considered in the time- 
evolution of the algorithm. Using the definitions of | i/j a ) 
and R a in Eqs. (22), (27), the inner product of the start¬ 


ing state and the perturbed eigenvectors can be expressed 


as 


(start| il) a ) 


1 / V3 

E a V | F' ( E a ) | 


{start | a 0 , {3 a ) , 


(40) 


where {start\ a 0 ,/3 0 ) is the overlap of the starting state 
with the marked vertex state. Applying our previous 
approximations for E + and F' (A + ), that is, for our per¬ 
turbed eigenstates closest to the Dirac point, we find 

| ( start\ ip±) | = + O • ( 41 ) 


Thus, our starting state is indeed a superposition of the 
perturbed eigenstates, \ij)±) as assumed in the previous 
section. This makes it possible to investigate the running 
time and success amplitude of the algorithm in more de¬ 
tail, that is, 


\{£\e~ im \start)\ 

VM-} — e iE+t (^| !>-)) 

= —~r |sin (E + t)\ . 

31 /| 


(42) 

(43) 

(44) 


where the perturbation matrices W s / t mark the source 
and target sites, respectively. 

As prescribed in Hz], the communication protocol op¬ 
erates by preparing a state localised on the source per¬ 
turbation, in our case on the neighbors of the source site, 
and allowing the system to evolve under our communi¬ 
cation Hamiltonian. The system then localises on the 
neighboring vertices of the target site. The effect of the 
two perturbation matrices W s / t amounts to disconnect¬ 
ing the two vertices from the underlying lattice. 

As noted before, our search algorithm has several dif¬ 
ferent optimal starting states depending on the position 
of the marked vertices. Therefore, we divide our commu¬ 
nication analysis into three different cases: the two sites 
i) are on the same sublattice and have the same optimal 
search starting state; ii) are on the same sublattice but do 
not share the same optimal search starting state; iii) are 
on different sublattices. Both of the first two cases can 
be treated by applying the reduced Hamiltonian method 
demonstrated earlier, but communication between differ¬ 
ent sublattices is not tractable by this method and so we 
focus on numerics in this case. 

For signal transfer between two sites on the same sub¬ 
lattice, we will assume, without loss of generality, that 
our communication takes place between two sites on the 
A-sublattice, (a s / t , P s /t) j here, the subscript s or t de¬ 
notes the source or target vertices respectively. Using ar¬ 
guments as employed in Sec. |IV| we can reduce the num¬ 
ber of relevant states. Ultimately, the search dynamics 
takes place in the subspace spanned by the basis | K) A , 

| A') , the A-type Dirac states, and |£ s ), |£ t ), the uni¬ 
form superpositions over the neighbors of the source and 
target vertices. This basis leads to the reduced Hamilto¬ 
nian 


It is clear from our earlier results for E± and I 2 that 
our algorithm localises on the neighbor state \l) in time 

T = = O (^y/NlnN^J with probability amplitude 

O ^1/vdn N^j. This confirms the logarithmic correction 
observed numerically and displayed in the inset of Fig. [3] 


VI. COMMUNICATION 

It has been demonstrated in HU that discrete-time 
search algorithms can be modified to create a commu¬ 
nication protocol. We show here that a communica¬ 
tion setup can be established also in the continuous-time 
search algorithm with minor changes due to the subtleties 
of our search. 


where 




( 0 

0 


eA 1 \ 

H = 


0 

0 

e iv ° 

e iVt 

V A 

e -i^s 

e~ iVs 

0 

0 




e -ivt 

0 

0 


= 


/t + 2 Ps/t) 

and 


, (46) 


L's/t — 


qf (2 ot s /t + Ps/t) ■ 

In the first case considered we assume that there 
are two perturbations on the graphene lattice located 
at the points ( a s ,/3 s ) A and (a t ,ft) A , chosen such that 
e*T( cl «+2^) = e *Tr(a t +2/3 t )_ This implies that a search 
for either vertex, using our search algorithm from 
Secs.[TV]&[V] would use the same optimal starting state. 
As the phases are equal, we drop the subscript on the per¬ 
turbation coordinates in what follows. Fixing our phases 
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FIG. 5. (Color online) Numerically calculated signal transfer 
on a 12 x 12 cell graphene lattice between equi vale nt vertices, 
using the communication Hamiltonian in Eq. (451. The sys¬ 
tem is initialised in | i a ) and localises on | It). Only the sum of 
probabilities to be found on the neighbor vertices is shown. 


and diagonalising the reduced Hamiltonian, we find it has 


eigenvalues Xf = ±2 ^ and Ag 2 = 0 with eigenvectors 

^±2) 

= \ (e* | K) a + e iv | K') A ± |4) ± 14)) 

(47) 


<*i) 

= ^(> m A -j v \£) A ) 

(48) 



= -L (14) - 14)) • 

(49) 


a single perturbation located at vertex (a s ,/3 s ) until it 
reaches maximum success probability, and then apply the 
second perturbation to the vertex (a t ,/3j)' 4 . The figure 
confirms the behaviour expected from the reduced model 
calculation. 

The communication mechanism essentially works in 
the same way as the quantum search algorithm, as it 
can be viewed as one marked vertex ‘finding’ another. 
The initial localised source state decays back towards the 
search starting state, and the system then searches for the 
target state. 

In our second case of signal transfer we analyse the be¬ 
haviour of a communication system where we have two 
perturbed vertices on the same sublattice, but with the 
restriction that e’T( a -+ 2 W -j- («t+2/3 t )_ As such, the 

two marked sites cannot interact via the same search 
starting state. However, the optimal search starting 
states are not orthogonal so that signal transfer is still 
possible, but the resulting interference effects make the 
analysis slightly more complicated. 

We rewrite the coordinates of the target in terms of 
the source, at = a s + x and /3* = /3 S + V- Again, fixing 
our pha ses and diagonalising the reduced Hamiltonian in 


Eq. (46), we find it has eigenvalues A^ = ±\/3y ^ and 


A±=± 


with eigenvectors 


Using these eigenstates, we may rewrite the source neigh¬ 
bor state as 


14) = 


1 

2 


(I*) - 


ip-: 



(50) 


Placing the system in the source state |4) and allowing 
the system to evolve under the reduced Hamiltonian, one 
finds 


|4) 

= l (e-^ |^ 2 ) 


t 


— l 

T' 





(51) 

i>-2 

))^ 

^0 

) 




(52) 


(«-/3) \K') 

A ) 


1 

+ 2 

(cos (A^t) 

- 

1)14 


(53) 


We note that, in the last line above, the term in the 
brackets involving only the Dirac states is actually the 
optimal search starting state for both perturbed vertices, 
defined in Eq. The system thus oscillates between 

the states localised on the neighbors of the perturbed 

vertices, |4) and |4), in a time T = f y^, via their 
optimal search starting state. 

Figure [5] shows the system evolving under the full com¬ 
munication Hamiltonian, Eq. ( [45] ) . The initial state used 
for the time evolution shown in Figure [5] is the true lo¬ 
calised state, that is, we run the quantum search with 


***) =S ( em ^ y> - 0 m 

oil's / \ 

+ ( e -^U+2y) _ 1 ) | k ') 

V’±i) = ~F~ 2 ~ (e i ~5 L(a:+2y) + l) I KY 


(54) 


1 -F \ 14) ± \ |4) 


(55) 


± 


e ii/ s e i^-(x+2y) 


\K!) A --\Q--\tt) . 


Using these eigenstates, one may write the source per¬ 
turbation as 


14) — 2 ( Vya) 


V’-Va)- 


V’l) - 


^-1)) ■ (56) 


The full expression for the time-evolution is rather cum¬ 
bersome; therefore, we only show the terms and prefac¬ 
tors we are interested in, namely |4) and |4) 


1^(4)=^ cos(A+ji) +cos(A+f) |4) (57) 


1 r 

- - [cos 


( A V3*) “ cos ( A ^) 14) ■ 


We can see here that the prefactors do not depend upon 
the coordinates of either the source or the target sites; 
the transport signal between sites on the same sublattice 
but with different optimal search starting state is thus 
independent of the position of the source or target site. 

In Figure [6] we show the system evolved under the full 
communication Hamiltonian. Again, the initial state is 
the true localised state on the nearest-neighbors of the 
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FIG. 6. (Color online) Upper: Numerically calculated sig¬ 
nal transfer on a 12 x 12 cell graphene lattice between non¬ 
equivalent vertices, using the communication Hamiltonian in 
Eq. ( |45| ). The system is initialised in \l g ) and localises on \lt). 
Only the sum of probabilities to be found on the neighbor ver¬ 
tices is shown. Lower: Analytically calculated behaviour for 
the same system, using the reduced Hamiltonian method and 
Eq. ([571. 


FIG. 7. (Color online) Numerically calculated signal transfer 
on a 12 x 12 cell graphene lattice between vertices on dif¬ 
ferent sublattices, using the communication Hamiltonian in 
Equation 45 The system is initialised in \l a ) and localises on 
| i t ). Only the sum of probabilities to be found on the neigh¬ 
bor vertices is shown. The two figures have the same source 
position but different positions of the target site. 


source vertex, obtained by running the search algorithm 
using one marked vertex until it reaches its peak success 
probability. The time-evolution is radically different to 
the previous communication case; the behaviour here is 
erratic with uneven peaks of probability at the two per¬ 
turbations involved in the protocol. However, there are 
still significant probability revivals. 

The transport behaviour from Eq. (571, calculated us¬ 
ing the reduced Hamiltonian, is also shown in Figure [6] 
The probability at time t = 0 has been scaled to match 
that shown in Figure [6] Our calculated behaviour has 
the same signal pattern as the numerically calculated 
behaviour from the full Hamiltonian, although over a 
shorter period of time. As our reduced model only makes 
use of the Dirac states and the perturber states, we lose 
the contribution to the time-evolution from the rest of 
the spectrum giving rise to logarithmic corrections as dis¬ 
cussed in Sec. 0 This leads to differences in the overall 
time scales for the reduced model and the full Hamilto¬ 
nian. This also supports our findings that signal transfer 
between all non-equivalent vertices on the same sublat¬ 
tice is the same. The communication protocol is again 
set up by a search mechanism in reverse, where one ver¬ 
tex finds another. The slightly erratic behaviour that 
emerges here is due to interference between the two sep¬ 
arate search mechanisms which interact due to the non¬ 
zero inner product of the three possible optimal search 
starting states for vertices on the H-sublattice. 

In our final case of signal transfer, we consider commu¬ 
nication between sites on different sublattices which we 
can not treat in the reduced Hamiltonian model. It has 


been demonstrated previously that perturbations to one 
sublattice do not interact with the Dirac states which live 
on the other sublattice. Therefore, when attempting to 
reduce the communication Hamiltonian as done before, 
we merely find that it decouples into two non-interacting 
reduced Hamiltonians, each describing a search proto¬ 
col on one sublattice. The interaction between the two 
sublattices is here facilitated due to interactions via the 
bulk of the spectrum. In what follows, we focus on nu¬ 
merics and inspect the behaviour for systems involving 
the same source perturber but different targets. The nu¬ 
merics show a finite number of different signal patterns, 
two examples are shown in Figure [7] We can see from 
these examples that the communication mechanism takes 
place over a much longer timescale with a superimposed 
oscillatory dynamics. 

The increase in timescale can be attributed to the weak 
nature of the interaction between the two perturbations 
due the fact that the localised states live mainly on one 
sublattice. The signal pattern in Figure[7]thus decouples 
into a fast oscillation between one site, say the source 
site, and the delocalised lattice state and a slow time 
scale on which a small amount of probability amplitude 
escapes into the other sublattice due to the weak interac¬ 
tion. This process continues until the recurrence proba¬ 
bility at the target perturbation reaches the same peak as 
the initial localised source state, and then the behaviour 
reverses. 
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VII. ALTERNATIVE METHODS OF MARKING 


In the previous sections, we have discussed a method 
of marking a vertex through modifying the hopping po¬ 
tential from the perturbed site to all three of its nearest- 
neighbors. Here, we will discuss alternative methods of 
marking a vertex which still keep the perturber interac¬ 
tion in the conic dispersion region of the spectrum. This 
can be done by altering the hopping potential in different 
ways as discussed above; this approach is necessarily a 
rank-2 perturbation and these two perturber states must 
meet at the Dirac energy. We will discuss several ways 
of applying such a perturbation in the next paragraph. 
Another approach is based on coupling extra sites to the 
lattice; this set-up will be discussed at the end of this 
section. 

Focussing on hopping potential perturbations, we may 
consider perturbing the bonds to any number of nearest- 
neighbors to any strength. We have seen previously that 
a symmetric three-bond perturbation successfully creates 
a search. Here, we demonstrate a search based on per¬ 
turbing the hopping potential from a given site to only 
one of its nearest-neighbors. That is, we use the same 
search Hamiltonian as in Eq. (13) with the perturbation 
matrix 


w = I a 0 ,(3 0 ) A (a 0 , 

Po\ B + \a 0 ,/3 0 ) B {a 0 , Po\ A , 

(58) 

and eigenstates 



I’V = ^ 

(| a 0 ,Po) A - |a 0 ,/? 0 ) B ) 

(59) 

|Hy - 71 

(jao, Po} A + |a 0 ,A>) S ) - 

(60) 


Note that it no longer makes sense to speak of a single 
marked vertex. Rather, our perturbation marks both 
vertices (a 0 , /3 0 ) A ^ B simultaneously; thus, it may be more 
accurate to view our perturbation matrix as marking a 
single cell of the lattice. 

For this type of perturbation, the avoided crossing used 
to generate search behaviour is not necessarily at 7 = 1, 
see Figure [HJ where the spectrum of our search Hamilto¬ 
nian as a function of 7 for a 12 x 12 cell torus is shown. 
As stated before, our single bond perturbation is also 
a rank-2 matrix, and so again there are two perturber 
states approaching the spectrum from the negative and 
positive regions of the spectrum. Inspecting the region 
around the Dirac energy, see inset in Figure[8j we see that 
avoided crossings are formed by four states in total, the 
two blue (solid) curves and the two red (dashed) curves. 
At 7 = | there is an exact crossing between the red and 
blue curves indicating that these states are orthogonal. 

We can make this avoided crossing picture clearer by 
breaking down the eigenstates of the search Hamiltonian 
in terms of the symmetries of the lattice. In particu¬ 
lar, we use the rotation operator C 2l which is a rotation 
of 7r about the mid-point of the perturbed bond. Con¬ 
sidering the action of C 2 on the marked vertices and 


y- 
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FIG. 8 . (Color online) Spectrum of single-bond search Hamil¬ 
tonian in Eq. (131 using perturbation from Eq. 58] as a func¬ 


tion of 7 for a 12 x 12 cell torus (N = 288). The spectrum 
is symmetric around tr> = 0. Inset: States nearest the Dirac 
energy have been coloured depending on their parity with re¬ 
spect to the C 2 operator: even (solid blue), odd (dashed red), 
undefined (dot-dashed black). 



FIG. 9. (Color online) Search on 12 x 12 cell graphene lat¬ 
tice with single-bond perturbation, using an optimal starting 
state. The behaviour at both marked vertices is the same, as 
is the behaviour at each of their neighboring sites. 


the eigenstates of the perturbation matrix, one finds 
C 2 \a 0 ,p 0 ) A/B = \a 0 ,p 0 ) B/A , C 2 \Wa) = -| W g ) and 
C 2 |W e ) = |We). The inset of Figure [8] shows the states 
nearest to the Dirac energy, where the same colour in¬ 
dicates the same C 2 parity. It now becomes clear that 
the states near the Dirac point actually form two avoided 
crossings, one for each parity with a minimum energy gap 
at 7 = |. 

Without going into details, we can again find optimal 
starting states using a reduced Hamiltonian approach 
(one for each of the avoided crossings); the evolution of 
the system using these optimal states is shown in Figure[9] 
for a 12 x 12 cell lattice. We find there is a significant 
localisation on the marked vertices and their nearest- 
neighbors, with the probability of being found on either 
of the marked vertices peaking at around 16 — 18%. We 
can also see that the probability of being found on each 
of the nearest-neighbors peaks at around 8%, resulting 
in a total probability of being found on the marked ver¬ 
tices and their nearest-neighbors of approximately 48%. 
The success probability fluctuates, as the probability am¬ 
plitude oscillates between the marked vertices and their 
nearest-neighbors. This is due to the probability ampli- 
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tude being constrained in the local area by the increased 
hopping potential between the two marked vertices. 

Although our demonstration is for a 12 x 12 cell lattice, 
numerical investigations show that the search behaviour 
remains the same as the lattice size increases with critical 
value fixed at 7 = |. We also find that the gap at the two 

avoided crossings scales as O (l/y/N In , in the same 
way as for the three-bond perturbation shown in Figure[3] 
As the search time is inversely proportional to the energy 
splitting at the avoided crossing, it also gives an estimate 
of the running time of the search T = O (^/N In N^j . 

We now turn to coupling additional sites to the lattice 
as a way of introducing a perturbation as also considered 
in [ 12 ] ; such a treatment is in many ways closer to exper¬ 
imental realisations as defects due to additional add-on 
atoms are quite common in graphene [23] . We focus here 
on the idealised case where an additional site is coupled 
to a single lattice vertex. We use the perturbation matrix 

w ( 7 ) = - I OLo, /3o) A (site I - | site) (a 0 , /3 D | A +7 I site) (site \ 

(61) 

for coupling the additional site | site) to the A-type vertex 
( a 0 ,(3 0 ) A and 7 is a free parameter related to the on¬ 
site energy of the additional vertex. The choice of the 
coupling terms is such that the binding energy between 
the additional site and the lattice vertex is the same as 
the internal couplings in the lattice. Thus, our search 
Hamiltonian is of the form 


H = -A + W ( 7 ). 


(62) 


In Figure [TUJ the spectrum of the search Hamiltonian 
is given as function of 7 . As our free parameter 7 only 
changes a single term, our spectrum only has a single per- 
turber state. One finds a clear avoided crossing around 
the Dirac energy (E = 0) when 7 = 0, that is, when the 
on-site energy of the additional site matches the on-site 
energy of the lattice vertices. 


In the reduced Hamiltonian picture, we write Eq. 62 
in terms of a basis consisting of the Dirac states and our 
perturber state, \site). Only three states are involved in 

the search, that is, {| K_) A , | K') , |szte)}, and H reduces, 
up to a different prefactor, to the same 3x3 matrix found 


in Eq. (141 for the three-bound perturbation. Thus, our 


previous reduced Hamiltonian analysis holds for this case 
and the optimal search starting states are the same; the 
time-evolution for the additional site search is shown in 
the inset in Figure [To] The only differences between this 
case and the three-bond marking are that the system 
localises on the additional site, and also the change in 
prefactor leads to a search time of T = jVN. 

We note that other, more realistic types of perturba¬ 
tions involving additional sites coupling to a lattice site 
and its neighbors, can be shown to result in effective 
search protocols. Also, the single additional site pertur¬ 
bation described here can, like the three-bond perturba¬ 
tion, be used to setup a communication protocol, with 



y 


FIG. 10. (Color online) Spectrum of the Hamiltonian in 
Eq. (621 as a function of 7 for a 12 x 12 cell torus (N = 288). 
The symmetry of the graphene spectrum is broken by the 
choice of the perturbation in terms of an additional site. Inset: 
Search on 12 x 12 cell graphene lattice with a single additional 
site perturbation, using starting state |s) from Eq. (17 1 


the same signal patterns and behaviours as found previ¬ 
ously. The initial state is in this case completely localised 
on the additional site. 


VIII. ALTERNATIVE NANOSTRUCTURES 

So far we have described the dynamics of searches on 
graphene lattices with periodic boundary conditions, that 
is, graphene on a torus. Here we consider more realis¬ 
tic boundary conditions such as nanotubes and graphene 
sheets. 

For the dynamics on nanotubes, we move to peri¬ 
odic boundary condition along one axis only and impose 
Dirichlet boundary conditions along the other directions. 
The properties of nanotubes are well-known [24], and it 
has been shown that the band structure of armchair nan¬ 
otubes, that is, nanotubes with armchair boundaries, al¬ 
ways allows for an energy at the Dirac energy regardless 
of the nanotube diameter; we will focus on these types of 
nanotubes. We are interested in searching on finite length 
nanotubes, where the band structure becomes discrete; 
the length of the nanotube are in addition chosen such 
that there exists an eigenenergy at the Dirac energy. An 
example of the cell we use to construct finite armchair 
nanotubes is shown in Figure El We choose the finite 
length of the nanotube to be along the horizontal axis 
and we close the underlying graphene lattice into a nan¬ 
otube along the vertical axis. 

We denote the basis states of sites in our nanotube as 
|to, A/B, l), where m indicates the m th A/B -type vertex 
in the horizontal direction in the I th cell. It is simple to 
see that, due to the periodic boundary conditions along 
the circumference of the tube, the vertical component of 
the eigenstates must be Bloch states and the horizontal 
component of the amplitudes must be sinusoidal in na¬ 
ture. That is, the eigenstates of the finite nanotube are 
standing waves along its length. 

By working through the tight-binding model for this 
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FIG. 11. Example of a armchair nanotube cell. The nanotube 
is periodic along the vertical axis and finite in the horizontal 
direction with a width of N x sites. 

system, one can show that (k x ,k y ) = (4^,0) is the only 
point where there exists an eigenenergy equal to the Dirac 
energy. It is also possible to show that the spectrum 
includes Dirac points when the number of sites along the 
length of the nanotube N x = 3r — 1, where r is an integer. 

As there is only one potential Dirac point for finite 
armchair nanotubes, it follows that there are only two 
Dirac states (one from the bonding and the anti-bonding 
regions of the spectrum). As we have k y = 0 at the 
Dirac point, the Bloch wave around the circumference of 
the nanotube is simply a uniform superposition. Another 
important feature of the Dirac states on the nanotube, 
which we have not encountered previously, is the exis¬ 
tence of nodal points where the amplitude of the eigen¬ 
state is 0. We find these nodal points occur at every third 
site along the horizontal axis. 

We focus on applying our three-bond perturbation 
from Sections EY1&IY] to the armchair nanotube. In our 
new labelling our perturbation takes the form 

W = I Too, A, l a ) {{m Q + l,B,lo\ + (m 0 - 1, B , l a \ 

+ (m 0 , B,l 0 \) + h.c.. (63) 

We have assumed that we are perturbing an A-type ver¬ 
tex on an even horizontal coordinate, so that the per¬ 
turbed site and its nearest-neighbors remain within one 
nanotube cell. While it is easy to re-write this per¬ 
turbation matrix for other sites we restrict ourselves 
to this form for simplicity. Note even in this form 
it can still be expressed in terms of a marked site, 
\m 0 , A,l a ), and a state which lives on the neighbors, 
V) = (l m o + 1, B, lo) + | Too - 1, B, lo) + | Too, B, l 0 )). 

Our search Hamiltonian for the nanotube is the same 
as in Eq. Inspecting the spectrum of the search 

Hamiltonian as a function of the free parameter 7 where 
the perturbation is not located on a nodal point of a 
Dirac state, one finds an avoided crossing around the 
Dirac energy when 7 = 1 . This spectrum is very similar 
to the one shown in Figure [3] for the search on the torus 
and will be omitted. If the perturbation is located on 
a nodal site the picture is different, however; one finds 
that there is an exact crossing at the Dirac energy when 




FIG. 12. (Color online) Searches on an armchair nanotube 
(N = 320) with triple-bond perturbation located in the bulk 
(upper figure) and placed near the edge of the nanotube (lower 
figure), using numerically calculated optimal starting state: 
\m 0 — l,B,l 0 ) (dotted blue), \m a + 1, B, l„) (dashed red), 

| m 0 ,B,l 0 ) (dot-dashed green), sum of neighbor probabilities 
(solid black). 


7 = 1 and searching is not possible. Recall that the effect 
of the perturbation at the critical point is to completely 
remove the site from the lattice, if the site already has 
zero amplitude then the perturbation will not interact 
with the Dirac state. 

Similar to the previous section, we numerically reduce 
the full search Hamiltonian in a basis consisting of the 
two Dirac states and the state living purely on the neigh¬ 
bors, \t/. The perturbed site, | m 0 ,A,l 0 ), is decoupled 
from the lattice. Through this process, we find two pos¬ 
sible starting states, one for each sublattice. The starting 
states are weighted superpositions of the Dirac states over 
the sublattice containing the marked vertex. The nodal 
points are excluded from this treatment. 

Figure [12] shows the system evolution for two searches 
using the numerically found optimal starting states. One 
search has the marked site located in the centre of the 
nanotube, the other has the marked site positioned near 
the edge of the nanotube. The dimensions of the nan¬ 
otube have been chosen so that there are 320 sites, com¬ 
parable to the searches shown in earlier sections. Com¬ 
paring these searches to the three-bond perturbation 
search on the torus in Figure [4j we see that there is 
a marked difference in behaviour and success probabil¬ 
ity induced by relaxing the periodic boundary conditions 
along one axis. 

Using a perturbation located near the edge, the lower 
image in Figure [12] we see a reduction in success proba¬ 
bility by a factor of 2-3 and strong fluctuations in peak 
height when compared to the search on the torus. Search¬ 
ing in the bulk, shown in the lower image of Figure [12] 
displays behaviour closer to searches on a graphene torus, 
but again with significant fluctuations at each peak. We 
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FIG. 13. Examples of two finite graphene sheets, with dimen¬ 
sions in terms of primitive cells (N x ,N y ) = (4,4). Along the 
vertical axis of both sheets are armchair edges. The horizon¬ 
tal boundaries are formed by bearded edges (left) and zigzag 
edges (right). 


FIG. 14. (Color online) Spectrum of triple-bond perturbation 
search Hamiltonian in Equation [13] as a function of 7 for a 
10 x 10 cell bearded graphene sheet (N = 200). We focus only 
on the relevant section of the spectrum. 


propose that this effect is due to the reflection of proba¬ 
bility amplitude from the edges of the nanotube. This is 
supported by the changes in the interference pattern in 
the signal as the perturbation is moved across the lattice. 
Numerics demonstrate that additional site perturbations 
and communication protocols can also be used on nan¬ 
otube lattices. 


We now move to working on graphene sheets, that is, 
removing the periodic boundary conditions along both 
axis. A detailed account of a specific version of this set¬ 
up together with experimental results has been presented 
in m for graphene sheets consisting of armchair bound¬ 
aries only. Note that such a configuration can not be 
achieved on rectangular sheets; the simplest configura¬ 
tion has the form of a parallelogram, see m for further 
details. The advantage of such a configuration is that 
the form of the boundary does not admit so-called ‘edge- 
states’. In the following, we will look at the influence of 
these edge-states in more detail by considering rectangu¬ 
lar graphene sheets, such as in Figure |13| In addition 
to armchair edges, here along the vertical axis, these 
sheets have boundaries formed by bearded edges (left) 
and zigzag edges (right) along the horizontal axis. These 
boundaries support edge states, i.e. states localised along 
these edges, with an eigenenergy close to the Dirac en¬ 
ergy [25] [26]. In the following, we will investigate, how 
the existence of these edge states influences the search. 

As in the case of the finite armchair nanotube, the 
imposition of Dirichlet boundary conditions at an edge 
generates sinusoidal eigenstates. As a result, there are no 
extended (bulk) eigenstates at the Dirac energy due to 
the inability to equate the quantised momenta with the 
necessary points in fc-space. We thus need to find other 
extended eigenstates in the sea of edge states near the 
Dirac point to undertake a search in this set-up. 

We mark sites using the triple-bond perturbation such 
as in Eq. (13). Throughout this section we choose the 
dimensions in terms of primitive cells of the graphene 
sheets, ( N x ,N y ) = (10,10); a bearded edge sheet thus 
consists of N = 200 sites and a zigzag sheet is formed of 


N = 218 sites. The spectra of the search Hamiltonians 


for the bearded lattice is shown in Figure 14 constrained 
to the energy region of interest. The spectrum for zigzag 
edges is very similar and is not show here. One can clearly 
see an avoided crossing around the Dirac energy at 7 = 1, 
the critical value for this type of perturbation. There are 
several states very close to the Dirac energy; these are the 
edge states which are all non-degenerate. Therefore, it is 
not possible for us to construct a superposition of degen¬ 
erate eigenstates which is optimal for searching. Rather, 
our initial starting state must be a single eigenstates of 
the unperturbed Hamiltonian H 0 = —A. 


Using the search Hamiltonian from Eq. (13) and fix¬ 


ing 7 = 1 , we proceed by allowing the system to evolve 
after being prepared in one of the unperturbed eigen¬ 
states. One finds for both types of lattice, that the edge 
states near the Dirac energy fail to produce localisation 
behaviour; the probability at the neighboring vertices of 
the marked site do not rise above noise levels. Rather, 
search behaviour only begins to emerge when we use the 
first delocalised, non-edge state as our initial state. We 
also note that we cannot search for sites near the zigzag 
or bearded edges, where the edge states exist. Only as 
we move further into the bulk of the lattice or along an 
armchair-type edge, the localisation behaviour returns. 
Note that the edge states can be viewed as a kind-of 
one-dimensional system and, as we saw from the scaling 
argument towards the end of Sec. [TT} one-dimensional sys¬ 
tems imply an energy spacing between successive energy 
levels of E n+ i — E n = O (IV -1 ). Thus, the perturber 
state interacts with many states in a dense part of the 
spectrum and the search fails. 


In Figure[l5]we show the search behaviour arising from 
marked sites placed at different positions on a bearded 
graphene lattice (zigzag lattice types display similar be¬ 
haviour). The results are similar to those found for the 
nanotube searches. There is a slight increase in variation 
of signal pattern with position, and an increase in success 
probability maxima as we move towards the centre of the 
lattice from either of the bearded edges. 
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FIG. 15. (Color online) Searching on a bearded edge graphene 
sheet with dimensions (TV X , N y ) = (10,10) using the Hamilto¬ 
nian, Eq. (131. All the searches are initialised in the first un¬ 
perturbed, non-edge eigenstate above the Dirac energy. The 
locations of the marked vertices are: upper - halfway along 
the left armchair edge, middle - centre of the sheet, lower - 
at the mid-point of the lower bearded edge and a third of the 
way along the vertical axis. (Colors and linestyles as in Figure 
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IX. DISCUSSION 


We have shown that continuous-time quantum search 
can be done effectively on a two-dimensional graphene 
lattice without the use of internal degrees of freedom. 
This is achieved by making use of the conical (lin¬ 
ear) dispersion relation in the graphene spectrum. The 

search succeeds in time T = O ^\/N In with prob¬ 

ability 0(1/ In TV). This is the same time complexity 
found in [31 for discrete-time searches and in [3] for 
continuous-time searches. To boost the probability to 
0(1), C(lnTV) repetitions are required giving a total 

time T = O N^j. Amplification methods [25]- 

[30] may be used to reduce the total search time further. 

Our main result focusses on perturbations which in¬ 
volve altering the hopping potential from a marked site 
to all three of its nearest-neighbors equally. We have also 
demonstrated other types searches based on perturbing 
the hopping potential in a cell and the adding extra sites. 
We have shown that search mechanisms can be utilised 
for the purposes of signal transfer. 


Our findings point towards applications in directed sig¬ 
nal transfer, state reconstruction, or sensitive switching. 
This opens up the possibility of a completely new type 
of electronic engineering using single atoms as building 
blocks of electronic devices. Our results demonstrate 
that a range of nanostructures constructed from graphene 
could be used to this end. 
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Appendix A: Calculation of Z 2 ( Sk/k l) 


We give here some details regarding the derivation of 
Eq. ( [35| and the logarithmic divergence of / 2 . It is clear 
that the dominant contributions to the / 2 fc summations 
come from the vicinity of the Dirac points. Approximat¬ 
ing the spectrum close to the Dirac points, one has 


I 2k =2V3N k 


-1 


E --- 

{p,q)&L {Sk,11P 2 + 25ft-,12P<7 + 5ft-,22? 2 ) 

1 0 ( 1 ) 


y - 

(p,q)eL (Sk^mP 2 + 2Sic,i2Pq + 5^ j22 g 2 )' 


(Al) 


Here the sums over integers p and q is over a rectangu¬ 
lar region L of the lattice Z 2 which is centered at (0, 0) 
and has side lengths proportional to y/N - the center 
(0, 0), corresponding to the relevant Dirac point, is omit¬ 
ted from the sum. As stated in the main text, for k > 1 


the corresponding sums converge which proves Eq. (341. 


For k = 1 we will establish constant C\ and C 2 such that 


Ci In TV < y 


1 


, , 5ft- iip 2 + 25ft- i2P<7 + 5ft- 22 q 2 

(p,q)£L — — — 


which then directly leads to Eq. (351. To establish Ci 

is positive 


< C 2 In TV 


(A2) 


note that because each term in the sum 
its value decreases by restricting it to a square region 
—ai'/N < p < aw/~N, — a\\fN < q < a±y/N which is 
completely contained in L. Up to an error of order one 
the sum over a square region can in turn be written as a 
sum over eight terms of the form 


fll y/~N p 

EE 


~i i 5ft: np 2 ± 25ft- \ 2 pq + Sk,2 2<? 2 

i=l 0=1 -’ - -’ 


(A3) 
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For fixed p we can find g m ax such that 
1 


E 


Sk,hP 2 ± 2Sx,i2Pq + Sx,22q‘ 
P 


> 


7max 


(A4) 


Sk,11P 2 2S'if,12Mmax + Sx,22q 2 
We may choose g max = b±p for some constant b i > 0, so 
ai\/iV p x/iV 


1 


^ ^ SWaiP 2 ± 2Sx,i2Pq + Sx,22q~ 


> c 


E 

P =i 


which diverges as In AC 


1 

P 

(A5) 


Establishing C 2 and the corresponding logarithmic 
bound from above, and thus proving Eq. (35), follows 


the same line by first extending the sum to a square of 
side length 2a 2 ^fN that completely contains L and then 
establishing 


E 


1 


Sk.up 2 ± 2Sx,i2pq + Sx,22q 2 
p 


< 


Sk,11P 2 ± 2Sk, 12Wmin + Sr ,229mm 


with q min = b 2 p. 


(A6) 
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